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Abstract 

We introduce mesoscopic and macroscopic model equations of chemotaxis with anomalous sub- 
diffusion for modelling chemically directed transport of biological organisms in changing chemical 
environments with diffusion hindered by traps or macro-molecular crowding. The mesoscopic mod- 
els are formulated using Continuous Time Random Walk master equations and the macroscopic 
models are formulated with fractional order differential equations. Different models are proposed 
depending on the timing of the chemotactic forcing. Generalizations of the models to include linear 
reaction dynamics are also derived. Finally a Monte Carlo method for simulating anomalous subdif- 
fusion with chemotaxis is introduced and simulation results are compared with numerical solutions 
of the model equations. The model equations developed here could be used to replace Keller-Segel 
type equations in biological systems with transport hindered by traps, macro-molecular crowding 
or other obstacles. 
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I. INTRODUCTION 



Diffusion and chemotaxis are fundamental to the motion of bacteria 111 , the directed mo 
tion of neutrophils in response to infection 2\ , hypoxia stimulated angiogenesis [3( and many 
other biological transport processes 2J. These transport processes can further be compli- 

n n 

cated by traps |4|, macromolecular crowding [5] or other obstacles resulting in anomalous 
subdiffusion characterized by an ensemble averaged mean square displacement of diffusing 
species, (r 2 (t)), that scales sublinearly in time, i.e., (r 2 (t)) ~ t 7 with < 7 < 1, p-[l7]. 
In this paper we introduce mesoscopic and macroscopic models for transport in biological 
systems with chemotaxis and anomalous subdiffusion. 

The classic macroscopic model for the evolution of a diffusing species, with concentration 
n(x, t), in the presence of a chemoattractant, with concentration c(x,t), is the Keller-Segel 
model Q 



dn 
~dt 



D 



d n 







dx 2 ^ dx 



dc 
dx 



n- 



(1) 



where D and x denote the diffusion coefficient and the chemotactic coefficient respectively. 
In this model if the chemoattractant is removed the evolution corresponds to standard 
Brownian diffusion with (r 2 (t)) ~ t. 



Anomalous subdiffusion can be modelled as 
or Continuous Time Random Walks (CTRWs) 



Pactional Brownian motion (fBm) 19M21] 



22 



sities 



23) with long-tailed waiting-time den- 



131 ] . Both of these models are non-Markovian and both exhibit the same sublinear 



scaling for the ensemble averaged mean square displacement. However the second moment 
of the velocity scales differently in the two models 24j and the time averaged mean square 
displacement differs from the ensemble averaged mean square displacements in the CTRW 
model, but not in the fBm model |25j. Both possibilities should be considered when inter- 
preting results from experiments using single particle tracking 25] or fluorescence recovery 
after photobleaching 26[ and a simple test has been devised for analysing experimental data 
to determine which model is most appropriate 27]. 

At the macroscopic level, anomalous subdiffusion can be modelled through a modified 
diffusion equation 



dC 
~dt 



P( 7 ,t)V 2 C 



(2) 



with the diffusion constant replaced by a fractional temporal operator. In the case of frac- 



2 



tional Brownian motion (fBm) this operator is given by 20|, |21] 



V I { 1 ,t) = D{ 1 ) 1 i 



7-1 



(3) 



n the Continuous Time Random Walk (CTRW) model 22j, with power law waiting times 



131 ] . the fractional temporal operator is given by 



Djj(7,f) = D( 7 ) 



at 1 " 7 



(4) 



where -D( 7 ) is a generalized diffusion coefficient with units of m 2 s 7 and 



51-7 



Y(t) 
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(5) 



at 1 - 7 y ' r( 7 ) at 7 (* - i') 1-7 

defines the Riemann-Liouville fractional derivative of order 1 — 7. There have been various 
attempts to modify the fractional macroscopic equations to include force fields and reactions 



The Fokker-Planck equation for diffusion in a force field can readily be generalized by 
replacing the diffusion coefficient with a time dependent fractional operator as above. This 
has been justified within the framework of CTRWs, for force fields that vary in space but not 



time 
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29] and for force fields that vary in time but not space 30] . However these deriva- 



tions do not extend to the more general case of anomalous subdiffusion in a general external 
force field f(x,t) that varies in both time and space. Two obvious possible generalizations 
in this case are 1311 



and 



32 



34] 



dn 
~dt 

dn 

at = 



DJJ 2 n - — ^-V (f(x, t)n(x, t)) 



dt 1 -^ 1 

at 1 - 7 7 



r/ 7 dt 1 -"! 



1 / a 1-7 



(6) 



(7) 



If the force field is purely space dependent then the two models are equivalent and the so- 
lution is time subordinated to the concentration of diffusing species in the standard Fokker- 
Planck equation. This temporal subordination is not physically appropriate for time de- 



pendent external force fields 



32]. However an alternate formulation using an Ito stochastic 



differential equation has been proposed with a modified subordination in which the force 



varies in real time rather than the random time 



331 ] . In chemotaxis there may be a physical 



link between the time scale of the diffusion and the time scale of the effective force field 



3 



since the latter depends on the concentration of another diffusing species. Similarly in the 



35 



36| the force field from the membrane 



fractional Nernst-Planck equation considered in 
potential depends on concentrations of the diffusing species. 

In Section II we introduce four different models of chemotaxis with anomalous subdiffu- 
sion. The different models are characterized by differences in the nature of the anomalous 
diffusion (fBm or power law CTRWs), and differences in the details of the underlying random 
walk processes. In Section III numerical solutions of the associated discrete space equations 
are obtained for each model. The numerical results are compared with Monte Carlo random 
walk simulations, with chemotactic forcing, on the same grid and using the same parameters. 
Differences between the model results are discussed in Section IV. 

II. FRACTIONAL CHEMOTAXIS DIFFUSION MODELS 

A. Model I 

To model chemotaxis with fractional Brownian motion we consider an ad hoc model in 
which we replace both the diffusion coefficient and the chemotactic coefficient by fractional 

temporal operators as in Eq. (j3J). This yields 

dn ^ d 2 n d ( dc\ 

dt 1 dx 2 ^ dx \ dx ) _ 

where 7 is the anomalous diffusion exponent, D 1 is the anomalous diffusion coefficient (with 

units m 2 s~ 7 ), and x-i is the analogous anomalous chemotaxis coefficient. This model equa- 
tion reduces to the standard Keller-Segel chemotaxis equation, Eq. (CQ), when 7 = 1. 

B. Model II 

A simple model for chemotaxis with fractional diffusion from power law CTRWs starts 
with the master equation 

t 

rii{t) = 7ii(0)$(*) + / (p r (zi_i,t')ni_i(*') +pi{xi +1 ,t)n i+1 {t)}^(t-t) dt (9) 







where ip(t) is a (power law) waiting time density, 

$(t) = I ^{t) dt (10) 



is the corresponding survival probability, and p r (x,t) and pi{x,t) are the probabilities of 
jumping from x to the adjacent grid point to the right and left directions respectively. 
These probabilities are dependent on the chemoattractant concentrations, c(x,t), at the 
neighbouring points of the point x at time t. The master equation, Eq. <M), is a continuous 



time representation o: 
Following Stevens 



the transition probability law in 37 1. 



37j, the probabilities of jumping to the left or right direction are based 



on the proportion of the chemoattractant on either side of the current point via 



and 



Pi(x i} t) 



Pr y^i ; t) 



v(x i+ i,t) 



(11) 



(12) 



v(xi-!,t) + t>(x m ,t)' 

where v(x, t) is a sensitivity function that depends on the concentration of the chemoattrac- 
tant: 

v(x, t) = exp (j3 c(x, t)). (13) 



Note that with the above we have 



Pl{x h t) +p r (x h t) = 1 



and 



Pi{x h t) -p r (x h t) 



Pc(Xi-!,t) _ J3c(x i+ i,t) 



(14) 



(15) 



e Pc(xi-i,t) _|_ e /3c(x i+1 ,t) ' 

Using the notation £{/(£)} (s) or f(s) to denote the Laplace transform with respect to 
time of a function f(t) we have the Laplace transform of Eq. 

ni(s) = 7i<(0)8(s) + {C{pr(xi-i,t)m-i(t)} (s) + C{ Pl (x i+1 ,t)n i+1 (t)} (s)}V>(a). (16) 



Using the identity 

S(a) = (l -?(«))/«, 
which follows from the Laplace transform of ( flOl) . we have 



(17) 



srii(s) - 71; (0) 



{-rii(s) + C{p r (xi-i,t)ni-i(t)} (s) + C {pi(x i+1 ,t)n i+1 (t)} (s)} 



(18) 
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We now consider a heavy-tailed waiting-time density which behaves for long-times as 

-1-7 



, / x Kit 

m ~ - - 

T \T 



(19) 



where 7 is the anomalous exponent, r is the characteristic waiting-time, and k is a di- 
mensionless constant. Using a Tauberian (Abelian) theorem [38[ we can write the Laplace 
transform for this density function as (for small s) 



^0) 



«r(i- 7 ) 

7 



ST 



(20) 



Using Eq. ( FlTl) . we then find the corresponding asymptotic form for the survival probability 

„ ^^-^ ^t-i (21) 
7 

and the ratio 



A, 



(22) 



where 



A, 



7 



«r(i- 7 )' 



Specific cases of waiting time densities are the Mittag-Leffler density 39 1 



d 
It 



T 



(23) 



where E 1 (z) is the Mittag-Leffler function 40], and the Pareto law used by 41| 

j/t 



[1 + t/r) 



1+7" 



(24) 



The corresponding values for A can be shown to be 



1 



(25) 



for ( 1231) and (1241 respectively. Note the ratio in Eq.( f22j) is only valid long-times for the 
Pareto density, Eq. ( )24|) . whilst it is exact for the Mittag-Leffler density for all times. In 
addition, if 7 = 1 we do not use Eq. f l24"|) but instead use Eq. (123]) . 
With Eq. d22J, Eq. (PJ now becomes 

sni(s)-ni(0) = — ^— - — {-ni(s) + £{p r (x i -. 1 ,t)n i -i(t)} (s) + £ {pi(x i+1 ,t)n i+1 (t)} (s)} . 



(26) 
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Noting that the Laplace Transform of a Riemann-Liouville fractional derivative of order a, 



where < a < 1, is given by 



40] 



dt a ~ l 



t=o 



(27) 



(28) 



we can invert the Laplace transforms in Eq. (l26|) to obtain 
dn- A d 1 ^ 1 

= ~^~faP^ { _7l i(*) + Vr{Xi-x,t)m-x{t) +pi(x i+ i,t)n i+1 (t)} 

where we have ignored the last term in Eq. (I2"7j) . Numerical solutions of this discrete space 
fractional differential equation for Model II are considered in Section III. 

The spatial continuum limit of Model II can be obtained in the usual way by setting 
Xi = x and Xj±i = x ± Ax and carrying out Taylor series expansions in x. Retaining terms 
to order (Ax) 2 and using the normalization pi(x,t) + p r (x,t) = 1 we first find that 

d 

- m(t) +p r (x i _ 1 ,t)n i _i(t) +Pi(xi+i,t)n i+ i(t) = Ax— (n(x,t)\pi(x,t) -p r (x,t)\) 

ox 



+ 



Ax 2 d 2 



-n(x,t). 



2 <9x 2 

This simplifies further after carrying out Taylor series expansions in Eq. f[To"|) , to arrive at 



(29) 



Pi(x,t) -p r (x,t) 



— tanh I p Ax— I . 



-I3Ax 



dc 
dx 



We can now combine the results in Eq.( |29l) and Eq.( |30l) with Eq. ( |28i) to obtain 

dn d 1 ^ 1 



dt dt 1 -"/ 



A 1 Ax 2 d 2 n l , A^(3Ax 2 d ( 'dc(x.t) , , 
(x, t) n — — v n(x, t) 



2t"< dx 2 



r 7 dx V dx 



+ 0(Ax 4 



and then taking the limit Ax — > and r — > 0, with 
and 



A^Ax 2 
2rT 



A 7 /3Ax 2 

7^ 



we have 



<9n 



dt l -i 



d 2 n(x,t) _ d_( dcjx, t) 
7 <9x 2 Xl dx\ dx H[X,) 



(30) 
(31) 

(32) 

(33) 
(34) 



Equation (|34|) provides a useful approximation for the space and time evolution of the con- 
centration of an anomalously diffusing species that is chemotactically attracted by another 
species. In the CTRW master equation, Eq.flH]) for this model the probabilities to jump left 
or right are determined at the start of the waiting times. We will consider another model 
in Section D where the probabilities to jump left or right are determined at the end of the 
waiting times. But first, in the next section, we consider an alternate formulation, based on 
a generalized master equation approach. 



C. Model III 



In this section we follow the generalised master equation approach of [30|, |42|, |43j , extended 
to take into account the effect of the chemoattractant. To begin we write the balance 
equation for the concentration of particles, n, at the site i 

^ = J?(t) - Jr(t), (35) 

where jf are the gain (+) and loss (-) fluxes at the site i. We also have the conservation 
equation for the arriving flux of particles at the point i given by the flux of particles either 
leaving the site i — 1 and jumping to the right or leaving the site i + 1 that move to the left: 

Jt{t)=Pr{xi-i i t)Jr 1 (t)+ Pl {x i+li t)J^ 1 (t) (36) 

where pi(x,t) and p r (x,t) are given in Eqs. ( TTT1) and ( TT2l) . We can combine Eq. (1361) and 
Eq. (1351) to obtain an evolution law for the concentration purely in terms of the loss flux, 
viz; 

-^= Pr ( Xi _ h t)Jr_ 1 (t)+p l (x i+l ,t)Jr +1 (t) - J~(t). (37) 
The loss flux at the site i is given by 

t 

j-(t) = ip(t)ni(0) + j V(t - t)J+(t) dt. (38) 
o 

The first term represents those particles that were either originally at % at t = and wait 
until time t when they leave. The second term represents particles that arrived at some 
earlier time t' and wait until time t to leave. Here ip(t) is the usual waiting-time density 



S 



used in Model II. We can combine Eq. ( 135]) and Eq. (1381) to obtain 



jr(t) = ij(t)n t (0)+ / i/>(t-t') 



— , /. dnAt) 

J - { f > + -1r 



(39) 



and then we can solve for the loss flux using Laplace transform methods. The Laplace 
transform of Eq.( l39l) with respect to time yields 



J-(s) = </>(sK(0) + i){s) Jr(s) + sn^s) - 12,(0) 



which simplifies further as 



jr(s) = &n t (s). 



(40) 



(41) 



Now using the approximation in Eq. ( |22|) for a heavy-tailed waiting-time density and invert 
ing the Laplace transform we have 

A 7 d^'n^t) 



■Kit) 



(42) 



t~i dt l ~~< 

We now substitute the expression for the loss flux, Eq.( l42|) back into the balance equation, 
Eq.f l37|) . to obtain 



driiit) A 7 ( d l 7 ni_i(t) d 1 J n i+1 (t) d 1 ^n^t) 

J p r {Xi_ U t) — — — + p l (x i+1 ,t)- 



(43) 



Numerical solutions of this discrete space fractional differential equation for Model III are 
considered in Section III. 

The continuous space representation of Eq.()43l) is found by setting x, = x and x,±\ = 
x ± Ax so that 

dn(x,t) Ay f d l ^'n(x — Ax, t) 



T' 



p r (x — Ax, t)- 



+ Pl ( X + A X ,t f" n t + _ AX ' t) - \ ■ (44) 



The continuum limit representation can then be found by carrying out Taylor series expan- 
sions about x, similar to the steps used to reduce (I28I) to (134")) . This results in the equation 



dn d 1 7 d 2 n(x,t) d ( dc(x,t) d 1 7 n(x,t) 



dt dt 1 -/' 



dx 2 



dx V dx 



dt 1 -"/ 



(45) 



Model II and Model III are similar to the fractional Fokker Planck equations, Eq. ([6]) and 
Eq.((7j) respectively with forcing from the chemotactic gradient 9c ^^ ■ 
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D. Model IV 

We now re-consider the master equation for the CTRW model but with the jump prob- 
abilities calculated after the particle has waited and immediately prior to jumping. The 
master equation in this case is given by 

t t 
n i (t)=n i (0)^(jt)+p r (x i -. 1 ,t) J m-xtfMt-t^dt' +pi(x i+u t) J n i+1 (t)i/j(t-t)dt . (46) 

o o 

It is convenient to introduce the auxiliary function 

niiit) = [ ni(t')ip(t - t) dt , (47) 
Jo 

which has the Laplace transform 

mi(s) = ni{s)$(s). (48) 

The Laplace transform of Eq. ( I46|) with respect to time can then be written as 

(sni(s) -rij(O)) S(s) = -$(s)rii(s) + C {p r (xj_i, ^m^t)} (s) + £ {pi(x i+1 ,t)m i+1 (t)} (s) , 

(49) 

and after the inverse Laplace transform, 

J -^$(t-t')dt =m i (t)+p r {xi- 1 ,t)m i - 1 (t)+pi(x i+1 ,t)m i+1 (t). (50) 

Proceeding to the continuum limit with Taylor series expansions about x, similar to the 
steps in Model II and Model III, we obtain 

/ •(« - 0^ * - - (^™M>) + 0(Ax% (51) 

o 

and then using the auxiliary function definition in Eq. (j4"7]) we find 

,.dn(x,t) , Ax 2 f d 2 n(x,t') 



f ^ , ,.dn(x,t) i Ax f d nix.t) , , /. ,/ 

J * (i " ' »T * * — j ~ 1 > dt 



_ Ax 2 04- I [ n (x, t'U(t - t) dt \ + 0(Ax 4 ). (52) 

ox \ ox J 
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Asymptotic expressions for the convolution integrals in Eq. (!52|) can be obtained by con- 
sidering asymptotic expansions in Laplace space and then inverting. Thus we now consider 
the terms 

t 



C J/«( t - ( ')*!^)«tt'U) = •(.)£ 



w . (53) 



„ , . d 2 n(x,t') , . /. , / I . . 7.. ,(9 2 n(x,s) . . 

£ < I —^-^(t -t)dtUa) = ^—^t 2 , (54) 





t 



C< / n(x,t')ip(t -t) dt } (s) = ${s)n(x,s). (55) 



For long times (small s) we have, with the use of Eqs. ( 1201) and (J21 



8(a) ~ fZ-^ + O^- 1 ), (56) 
?(s) ~ l-^j^ + 0(s 7 ). (57) 

vi 7 

Using these expansions in Eqs. (I53p - (j55p and taking the inverse Laplace transforms to 
replace the convolution integrals in Eq. ( 1521) we obtain 

r7 d 7 " 1 ggfoj) ~ fr^^MM) _ Ax2/3 A / ^M n f x ^ + o(Ax 4 ) (58) 



and in the limit Ax — > and r — > 



dn(x,t) d 1 7 



£) d 2 n(M) 9 ( dc ^ n (x,t) 



(59) 



as previously in Eq. ( 1341) . 

Conversely, for short times (large s) we have 

S(s) ~ -+0(s-^- 1 ), (60) 
~ ^-^ + 0( S - 2 ^), (61) 

where z/ 7 = 7 and £> 7 = 1 if use the Mittag-Leffler ( 1231) and z/ 7 = 1 and £? 7 = 7 if we use the 
Pareto (1241) density. The resulting equation for Eq. (1521) becomes for short-times 



dn(x, t) _ a 1 ^ d 2 n(x, t) J 9 / dcpM) a v m(x,t) ' 



dt 7 dt l ~ u -< dx 2 ^ dx dt \ dx dt 
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with the modified coefficients 



and 



B At 2 

= (63) 

7 2t u i 



*7 = Vt " (M) 

In the case of the Mittag-Leffier density (1231) . the short-time equation can be simplified to 



dn(x,t) d 1 7 d 2 n(x,t) d fdc(x,t)d 1 J n(x,t) 

7 »*1 --v ^7 



at 7 at 1 - 7 ax 2 ^ 7 ax v ax at 1 - 7 



, a / a 2 c(x,t)a- 7 n(x,t) > .. r 
* 7 ax v axat at- 7 J ' 1 j 



Equation (163)1 is similar to Model Ill's governing equation given by Eq. ( 145]) if we note that 
the last term which will be close to zero near t = 0. Conversely if we use the Pareto density 
( 1231) then Eq. ()62]) becomes 



f ffc(M)a-MM) ' 
7X1 ax V a^at at- 1 J 1 J 



dn(x.t) ^ d 2 n(x,t) d fdc(x,t) 

"7^i^o 7X177- ^r—n(x,t) 



at a^ 2 ax v ax 



which is similar to the standard chemotaxis equation ([T]) (if we again consider the last term 
to be small) except for the presence of the 7 term. This term will, in effect, slow the initial 
temporal behaviour of the solution (linear rescaling). 

Note that if we use the Mittag-Leffier density in ( J23|) then we see that the mesoscopic 
equation Eq. f H6|) bridges the gap between Model II and Model III. At short times it recovers 
Model III, Eq.( T45|) . whereas at long times it recovers Model II The latter shows that for 
long times, compared to the characteristic time r, there is no difference between evaluating 
the chemotactic probabilities at the time before or after the particle waits. 

III. FRACTIONAL CHEMOTAXIS REACTION-DIFFUSION MODELS 

In this section we consider extensions of the CTRW based fractional chemotaxis diffusion 
models to incorporate reactions. In the absence of chemotaxis, extensions of CTRW based 



fractional diffusion models to include linear reactions were derived in 
to include nonlinear reactions were derived in [46], k' 



44 



45] and extensions 
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A. Model II 



Following the approach in 44) we can incorporate reactions in CTRW models by in- 
creasing or decreasing the concentration of particles during the waiting times by an amount 
proportional to the evolution operator for the reaction dynamics. The master equation for 
Model II with linear reaction dynamics incorporated in this way becomes 



n i (t) = e kt n i (0)$(t) + J \p r {xi-i, t'K-i(t') + Pi(x i+1 ,t')n i+1 (t')j eH* ^^(t-t^dt (67) 

o 

where k is the per capita rate gain (k > 0) or loss (k < 0) of particles. Again Laplace trans- 
form methods can be used to convert the integral equation representation into a (fractional) 
differential equation. The Laplace transform of Eq. ( 167|) with respect to time yields 



Ui(s) = 7^(0) $(s -k)+C {PriXi^fjUi-^t)} (s) ip(s - k) 

+ £ {pi(x i+1 , t)m+i(t)} (s) $(s - k), (68) 

and after rearranging we find 
sni(s) - n»(0) = krii(s) + — — {-n»(s) + C {p r (x i -. 1 ,t)n i -. 1 (t)} (s) 



$(s - k) 



+£{p l {x i+1 ,t)n i+1 {t)}{s)}, (69) 



where we have used (IT7|) . With the result in Eq. (l22|) we can inverting the Laplace transform 
to obtain 



drii 
~dT 



kt ^7 



-— (e kt \-rnit) +p r {x i - l ,t)n i „ l {t) + Pl {x i+l ,t)n i+l {t)}) + km{t) (70) 



where the Riemann-Liouville fractional derivative has been replaced by a modified fractional 



derivative 



45|. 



The continuum limit, found by taking Taylor series expansions about x, is 



dn 
~dt 



-kt 



d 2 n(x,t) d fdc(x,t) 
7 dx 2 ^ 7 dx \ dx n X) 



+ kn(x,t). (71) 



We note if n(x,t) is not self-chemotactic then the solution of ( I7ip is given by n(x,t) 
e kt y(x,t) where y(x,t) is the solution of (|34|) with y replacing n. 
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B. Model III 



To incorporate reactions in Model III we start by modifying Eq. ( 135]) to 

drii(t) 



dt 



J?(t)-Jr(t) + km(t) 



(72) 



where k, again, is the per capita rate gain (k > 0) or loss (k < 0) of particles. We also 
modify the expression for the loss flux, J^(t), in ( 1381) to 



J-(t) = e fc V(*H(0) + J eH*"* )^{t - t')J?(t') dt 

o 



(73) 



where the exponential factors take into account the per capita addition or removal of particles 
as in 44]. 

Now solving for the gain flux, J^it), in Eq. (!72|) we find 



JHt) = Jr(t) + ^-kn t (t) 



(74) 



and using Eq. ( 1731) gives 



jr(t) = e fe V(*K(0) + / e H'-* - f) 



Now using Laplace transform theory we find 



dt'. (75) 



Jj (s) = ij)(s — k)rii(0) + ip(s — k) J i (s) + srii(s) — rij(0) — krii(s) 



and upon solving for the flux, we find a similar expression to Eq. (T4T 



(76) 



J i ( s ) = 7T n i( s )- 



(77) 



$(s - k) 

This Laplace trasnform can now be inverted to find the loss flux given by the modified 
fractional derivative 44j, |45| of the concentration at i, 

,Av d 1 ^ 



■K{t) = e k ^—-{e- k %{t)) 



(78) 



where is a constant given by 1 or 1/T (1 — 7) if we use the waiting-time density in Eq. ( 123|) 
or Eq. ( |24|) respectively. 
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Now using Q3BJ and (178"]) in ([72]) we find 

-e kt ^^(e- kt Mt)) + k ni (t). (79) 

The continuum limit following from setting ± Ax, and Taylor series 

expansions about x, is given by 

dt~ V " e dt^{ 6 dx 2 J X "dx{ dx 6 dt^ [e H*,t)))+Kn{x,t). 

(80) 

C. Model IV 

The master equation for Model IV, Eq.f HoT) . modified to include linear reaction dynamics 
is given by 

t 

m{t) = e kt ni (0)^(t) +p r ( Xi - U t) [ n i -i{t')e k (*~ t ')il>(t-t')dt' 



+ Pl (x i+1 ,t) J n i+l (t)e k V 'J^t-t^dt. (81) 
'0 

Following similar steps used to simplify Eq. ( H6j) and using Taylor series expansions we find 

t ( t , 

k(t-t)^, ,.dn(x,t) ,1 Ax 2 f d 2 n(x,t) k(t~t) ,, \ , > 

A* l h{t-t dt ~ — J £j U K* )^(t ~t)dt 


+k j e k ^~ t ')^{t-t )n(x,t )dt -Ax 2 f5^- / n(x, - t) dt j +0(Ax 4 ). 

V / 

(82) 

Using Laplace transforms and the asymptotic expressions in Eqs. ([55]) and ([57]) (evaluated 
for s — k small) we arrive at Eq. ( 17 ip for long times. 

For short times we find 
dn(x,t) _ kt d 1 ^ ( _ kt d 2 n(x,t)\ kt d d ( dc(x,t) d~ v '< , _ kt *\ . . 

(83) 

with z/ 7 and the modified coefficients as defined previously for Model IV in section [III 
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IV. NUMERICAL SOLUTIONS AND MONTE CARLO SIMULATIONS 



It is straightforward to obtain numerical solutions of the above model equations using dif- 
ference approximations. In this section we describe numerical solutions for self-chemotactic 
variants of the model equations and we compare the solutions with Monte Carlo simulations. 

Implementation details for the Monte Carlo simulations are described in the Appendix lAl 
In the results reported here simulations were conducted on a one-dimensional lattice using 
the Pareto waiting-time density (1241) with the characteristic waiting-time r = 0.1, fractional 
exponent 7 = 0.5, and chemotactic sensitivity, (3. The simulation results are from an average 
of 200 runs with 10, 000 particles initially located at the origin. 

As our starting point for numerical solutions of the macroscopic models we consider the 
discrete space variants of Model II, III, and IV given by Eq f|28|) . fj43|) . and (T46]) respectively. 
A discrete space variant for Model I, analogous to the discrete space variant for Model II, 
is given by 

dnAt) A 7t 7 ~ 1 

—7— = — — - — \p r (xi-i,Xi,t)ni-.i(t) + pi(x i+1 ,Xi,t)n i+ i(t) - m(t)] (84) 
at t 1 

The discrete space equations for Models II and III were solved using an implicit time 
stepping method with the fractional derivatives approximated using the LI scheme 48] as 
in 49(. For Model IV, the integrals in the discrete space representation were approximated 
by taking the unknown concentration, rij±i(t), to be piecewise linear in time. 

The numerical solutions of the discrete space equations and the Monte Carlo simulations 
were performed using the same space grid size and similar values for the parameters r, 7, 
and (3. We also chose the same initial condition (one at the origin and zero elsewhere). The 
constant A 7 was chosen as in Eq. ( |25l) since we used the Pareto density, Eq. ( l24l) . in the 
simulations. 

In Figure [1] we compare the Monte Carlo simulation results with the numerical solution 
of the discrete space equations for each model with the chemotactic sensitivity parameter, 
(3 = 0.1. Further results are shown in Figures [2] and [3] for the sensitivity parameter values 
(3 = 1 and (3 = 10 respectively. 

The numerical solutions for Model III Eq ( H3l) and Model IV ( flol) are in close agreement 
with the Monte Carlo simulations at all times. The numerical solution for Model II does not 
fit with the Monte Carlo simulations well at short times but it provides a good fit at long 
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times (t = 20). The numerical solution for Model I does not fit the Monte Carlo simulations 
well, especially for small values of /3, where the predicted shape near the origin is smoother 
than that exhibited by the simulation data. 

The closer agreement between the numerical results for Models III and IV and the Monte 
Carlo simulations is due to the timing of the chemotactic forcing. In models III and IV, 
and in the Monte Carlo simulations the chemotactically influenced jumping probabilities are 
determined at the end of the waiting times, whereas in Model II they are determined at the 
start of the waiting times. This difference is less marked if the chemotactic concentration 
varies slowly in time. 

Overall, the numerical solutions for Model III provide better agreement with the Monte 
Carlo simulations than the numerical solutions for Model IV. This better agreement can be 
seen at the intermediate value of /3 — 1 in Figure |2j This better agreement may be due to 
differences in numerical errors in approximating the discrete space equations for Model II 
and Model IV rather than due to differences between the equations themselves. 



V. SUMMARY AND DISCUSSION 



The correct form of the fractional Fokker-Planck for particles undergoing anomalous 
subdiffusion in an external space and time varying force field has been an open problem. 
In the absence of a force field, subdiffusion can be modelled with a fractional temporal 
derivative operating on the spatial Laplacian. For subdiffusion in a purely space dependent 
force field the fractional temporal derivative can be put to the left of the standard terms 
on the right hand side of the standard Fokker-Planck equation 13j, |29|, |3l|, |50j as in Eq.(j6]). 



However the consensus has been that for subdiffusion in an external space-time-dependent 



force field the fractional temporal derivative should not operate on the force field [32N34j . 
as in Eq.fJH). The modelling is further complicated if the force itself is affected directly or 
indirectly by the subdiffusing particles. This is the case in fractional electro-diffusion and 
fractional chemotaxis diffusion, the case considered here. 

In this article we have introduced and investigated four models for particles undergoing 
anomalous subdiffusion in the presence of chemotactic forcing. The first being based on an 
adhoc extension to the fractional Brownian motion equation (Model I), two models based on 
Continuous Time Random Walk (CTRW) master equations where concentration-dependent 
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FIG. 1: One-dimensional sub diffusive chemotaxis simulation results (circles) with f3 = 0.1 showing 
the concentration profile at t = 0.4 (blue), t = 2 (red), t = 4 (black), and t = 20 (brown) compared 
with the numerical solution of the governing equations for models I (top left), II (top right), 
III (bottom left) and IV (bottom right) (solid lines). 
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FIG. 2: One-dimensional sub diffusive chemotaxis simulation results (circles) with f3 = 1.0 showing 
the concentration profile at t = 0.4 (blue), t = 2 (red), t = 4 (black), and t = 20 (brown) compared 
with the numerical solution of the governing equations for models I (top left), II (top right), 
III (bottom left) and IV (bottom right) (solid lines). 
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FIG. 3: One-dimensional subdiffusive chemotaxis simulation results (circles) with (3 = 10.0 showing 
the concentration profile at t = 0.4 (blue), t = 2 (red), t = 4 (black), and t = 20 (brown) compared 
with the numerical solution of the governing equations for models I (top left), II (top right), 
III (bottom left) and IV (bottom right) (solid lines). 
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jump probabilities were evaluated before (Model II) or after (Model IV) the particle waiting, 
and a fourth model derived from a generalized master equation (Model III) . Concentration- 
dependent jump probabilities were used to incorporate the effect of chemotaxis in discrete 
space representations of the models and in the Monte Carlo (MC) simulations. 

Evaluating the jump probabilities prior to waiting in the CTRW formulation (Model 

II) resulted in a macroscopic equation (valid in the long time limit) with the fractional 
derivative acting upon the chemotactic gradient. Conversely, using a generalized master 
equation approach with the probabilities evaluated after waiting but prior to jumping gave 
a macroscopic equation where the fractional derivative does not act upon the gradient (Model 

III) . The CTRW formulation with the jump probabilities evaluated after waiting (Model IV) 
could only be reduced to a Fractional Fokker-Planck continuum equation in the asymptotic 
limit for long and short times. For long-times Models II and IV coincide whilst for short- 
times we found Models III and IV coincide asymptotically if a Mittag-Leffler density is 
used. 

We also introduced Monte Carlo methods for simulating anomalous sub-diffusion in a 
chemotactic force field. In the Monte Carlo simulations the chemotactically influenced jump 
lengths were computed at the end of the waiting times, similar to Models III and IV. This 
could explain the excellent agreement we found between numerical solutions for Models III 
and IV and the Monte Carlo simulations. The numerical solutions for Model II also showed 
good agreement at long times. The numerical solutions based on the fractional Brownian 
motion equation, did not agree well with the Monte Carlo results. 

The fractional chemotaxis diffusion models were further generalized to incorporate linear 



reaction dynamics. As in previous research, 44], |45|, |5l| we found that the incorporation of 



linear reactions required the replacement of the Riemann Liouville fractional derivative with 
a modified version, in addition to including the linear reaction term. 

The fractional chemotaxis diffusion equations developed in this paper provide a new 
class of models for biological transport influenced by chemotactic forcing, macro-molecular 
crowding and traps. We have recently generalized these models to include arbitrary space- 



and-time dependent forces 



52|. 
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Appendix A: Monte Carlo Simulations 

In this section we briefly describe the Monte Carlo method used to simulate chemotaxis 
on a periodic one-dimensional lattice with long-tailed waiting-time density (subdiffusion). 
For each simulation run, the following steps are conducted 

1. Set the number of grid points, simulation time, and the initial number of particles. 

2. Initialise the parameters for the waiting-time and jump-length probability density 
functions. 

3. Set up the initial particle positions. 

4. For each particle generate a random waiting-time, St, (time of the first jump). 

5. Initialise the output time t out = At. 

6. Find the time of the next jump by finding the minimum of all jumping times, tj ump . 

7. If tj ump > t out then go to step ([H]) otherwise go to step Q. 

8. Store the current particle positions. Add At to t out . If t out exceeds the simulation time 
then simulation ends otherwise go to step ([7]). 

9. Generate a random jump-length, Sx (see below). 

10. Generate a new waiting time and update this particle's jumping time and position 

(tjump tj um p -\- St X new %old ~\~ Sx ) . 

11. Go to step © 
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1. Generation of Waiting Times 



The waiting times for each particle/jumper were generated by comparing a uniform ran- 
dom number, r G (0, 1), with the cumulative probability function of the waiting-time density. 
We use the Pareto density ff2fl§ as the density which has used by Yuste, Acedo, and Lindeberg 



4l|. The generated waiting-time is given as 

St = t ((1 -r)~T - 1) (Al) 



where r G (0, 1) is a uniform random number. 
2. Generation of Jump Distances 

The jump distance for each particle/jumper is generated by comparing an uniform random 
number, r G (0, 1), with the cumulative probability function of the jump- length density. 

For the simulations we use the jump-length probability density nearest neighbour jumps 
only: 

I -Ax, < r <p h 
8x = < (A2) 
i Ax, pi < r < 1 

where Ax is the grid spacing and pi = pi(xi, t) is the probability of jumping to the left given 
previously in Eqs. f lTTj) and f fT3]) . 

To evaluate the probabilities of jumping to the left or right for Eq. (1A2I) . requires the 
approximation of the chemoattractant concentration, c(xi,t), in Eq. (TH?]) . This is estimated 
by the proportion of chemoattractant particles at the grid point, x iy compared with the total 
number of particles in the system. 
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